Stochastic population dynamics in turbulent fields 



M. H. Vainstein 

Instituto de Fisica and Centra Internacional de Fisica da Materia Condensada, 
Umversidade de Brasilia, CP 04513, 70919-970, Brasilia-DF, Brazil 

J. M. RubS 

Departament de Fisica Fonamental, Universitat de Barcelona, Diagonal 647, 08028 Barcelona, Spain 

J. M. G. Vilar 

Integrative Biological Modeling Laboratory, Computational Biology Program, 
Memorial Sloan- Kettering Cancer Center, 1275 York Avenue, Box 460, New York, NY 10021, USA 

The behavior of interacting populations typically displays irregular temporal and spatial patterns 
that are difficult to reconcile with an underlying deterministic dynamics. A classical example is the 
heterogeneous distribution of plankton communities, which has been observed to be patchy over a 
wide range of spatial and temporal scales. Here, we use plankton communities as prototype systems 
to present theoretical approaches for the analysis of the combined effects of turbulent advection 
and stochastic growth in the spatiotemporal dynamics of the population. Incorporation of these 
two factors into mathematical models brings an extra level of realism to the description and leads 
to better agreement with experimental data than that of previously proposed models based on 
reaction-diffusion equations. 



I. INTRODUCTION 



Plankton patchiness has many causes, with origins both in biological and physical factors. Over small scales (1 mm 
to 10 m), biological factors on the individual scale such as mating, predator avoidance, finding food and the diel 
vertical migration (upward and downward swimming at certain times during the 24 h day) are crucial factors for the 
emergence of plankton patchiness On larger length scales (10 m to 100 km), physical processes such as turbulence, 
currents and eddies are the principal causes of patterns. 

It has been shown that the relative intensity of zooplankton patchiness is greater than that of phytoplankton 
at all spatial scales 0- Mackas and Boyd developed a method to count individual particles in a continuous 
stream of seawater. This permitted them to use shipborne sampling to find transects of zooplankton abundance 
and chlorophyll fluorescence, which provides an appropriate method for estimating phytoplankton concentration, and 
to make a spectral analysis of the spatial heterogeneity. They found that small spatial scale contribution to the 
total patchiness is much greater for zooplankton, which is reflected in the fact that power spectra are less steeply 
sloped for zooplankton than for phytoplankton. Moreover, they found that the spatial patterns of phytoplankton and 
zooplankton are negatively correlated. 

Satellite imagery allows the obtention of greater spatial coverage and a better time resolution than shipborne 
sampling. The analysis of changes in the spectrum of visible light due to absorption and fluorescence of chlorophyll 
pigments allow the estimation of phytoplankton biomass. Cower et al. 4] performed a spectral analysis of satellite 
images, the results of which lead them to believe that phytoplankton patchiness is controlled by mesoscale (10-100 km) 
water motions. Phytoplankton have small mobility, consisting of vertical migration limited to over a few meters per 
day and may therefore be considered almost as a passive scalar advected by the ocean currents. 

Many model systems have been put forward to explain pattern formation in such planktonic systems. To cite a few, 
Abraham Q has considered the role of non-diffusive advection in plankton pattern formation and Lopez et al. Q and 
Hernandez-Garcia et al. Q offer a comprehensive overview of modeling of plankton dynamics as chaotic advection. 
However, the important effects of noise have usually been neglected. 

This work is organized as follows: in the next section, we briefly discuss models of population growth used in 
ecology; in Sect. IIIII we describe the effects of noise in such models; in Sect. IIVI we analyze how the inclusion of 
diffusion and advection can generate patterns in population dynamics; in Sect.|Vl we discuss spectral analysis, one of 
the main tools used in studying plankton patterns; we give final remarks in Sect. I VII 
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FIG. 1: Image of a phytoplankton bloom in the Barents Sea north of Russia, captured by the Moderate Resolution Imaging 
Spectroradiometer (MODIS) on NASA's Aqua satellite on August 29, 2006. Source: NASA's Earth Observatory isj]. 

II. POPULATION GROWTH MODELS 

In ecology, continuous growth population dynamics is usually modelled by a system of differential equations which 
correspond to the conservation equations for each population, in which the rate of change of the different species is 
related to birth, death and migration processes. Predator-prey models are usually defined through 

— = F^(iV,P), and (1) 

^ = FpiN,P), (2) 

where N{t) represents the prey population, P{t) the predator population and the functions Fn{N,P) and Fp{N,P) 
account for the usually nonlinear interactions between the species. A general starting point is the Lotka-Volterra 
model, in which Fn{N, P) = N{a — bP), and Fp{N, P) = P{cN — d), with a, b, c and d positive constants. However, 
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this model has some drawbacks, since its sohitions are not structurally stable [9|, meaning that small perturbations can 
exert large effects on the amplitude of oscillation, for example. An improvement is given by a logistic growth of the prey, 
which limits the prey population even in the absence of predators. This feature is taken into account by the function 
Fn{N,P) = rN{l — N/K) — cPf{N), where r is a linear birth rate, K is the carrying capacity of the environment, 
and c a parameter related to the predation rate. The response function f{N), to be more realistic, should saturate 
to account for satiation. A possible choice is the HoUing type III functional response f{N) = 7VV(l + iV2) 0. An 
equation with these terms and considering a constant number of predators P{t) = P was used to model the outbreak 
of the spruce budworm (T3|. In that case, the carrying capacity is related to the foliage (food) available in the trees, 
and the functional response models a switching on of the predation rate by birds at a given threshold population of 
the prey. In other words, if the prey population is low, the predator (birds) will find a different source of food and 
if the prey population is high enough, it will become a source of food. Moreover, the saturation in the functional 
response represents the fact that there is a maximum uptake of food by the birds. 

A model similar to the one for the spruce budworm was used together with a time-dependent growth rate r = r(t), 
and Fp{N,P) = P{gf{N) — e), with g and e positive constants, to model the interaction of phytoplankton and 
zooplankton with different timescales that lead to rapid increases of phytoplankton population known as "spring 
blooms" and "red tides" p^ . 



III. THE ROLE OF NOISE 




FIG. 2: Population dynamics given by Eqs. Q and ((5j with r = 0.3, K = 4.0, c = 2.0, g = 0.1, e = 0.05. The populations 
fluctuate around the stable equilibrium Neq — 1.0 and Peg ~ 0.225. For larger values of the noise, the system can be driven 
to the unstable equilibrium in which the zooplankton population becomes extinct, (a) Phytoplankton concentration, (b) 
Zooplankton concentration. 

Deterministic population models rarely capture all of the features of dynamical systems, since real systems are 
influenced by changes in the environment which are in essence random. This inherent unpredictability of the envi- 
ronment may manifest itself in fluctuating birth rates, carrying capacities and other parameters which characterize 
biological systems. A deterministic treatment is more adequate for very large populations. However, the population 
of zooplankton is considerably smaller than that of phytoplankton and a stochastic approach is more realistic. The 
reason for the inclusion of noise is the fact that zooplankton interact with fish and whales which are present in even 
smaller numbers and are far from being evenly distributed. 

Nevertheless, in many cases, the inclusion of noise may destabilize a system and lead to the extinction of a given 
species. The linear stability of a system of m species 

^^F,{N,{t),...,NUt)) (3) 

depends on the eigenvalues of the matrix with components aij = {dFi/dNj)eq which governs the dynamics near 
equilibrium in the linear approximation. For a system without noise to be stable, it is necessary that all eigenvalues 
have negative real parts. If we define A as minus the largest real part of the eigenvalues, then the stability criterion 
becomes A > 0. In the presence of fluctuations characterized by variance cr^, the stability criterion should be changed 
to A > cr^ . The noise will generally tend to decrease the average population number by increasing the severity of 
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the population fluctuations. To see the effects of noise, we will consider the model 

=rN{\- NlK)-cPN^/{\ + N^), (4) 

— ^P{gNy{l + N')-e)+m, (5) 

where ^(t) is a random term, and N and P are the phytoplankton and zooplankton populations, respectively . As 
pointed out in , the intrinsic rate of zooplankton growth is lower than that of phytoplankton, what should lead to 
the formation and maintenance of only very large scale features. They concluded that a mechanism (maybe behavioral 
in origin) with shorter time scale should be at play to account for the small-scale patchiness observed. This is the 
motivation for the inclusion of a multiplicative noise term, since fluctuations originate from processes which depend on 
the local concentration, such as the reproduction process and the consumption by larger animals. Therefore, we take 
{^{t)£,{t')) = 2[crP{t)]'^S{t — t'), where <jP{t) is a fluctuating growth rate. An a posteriori reason for not considering 
additive noise is that it may allow the concentration of zooplankton to reach negative values, which is unreasonable. In 
the study of stochastic population dynamics in the Poisson approximation, an expression for the fluctuations around 
the deterministic limit in which the noise amplitude is state dependent is also arrived at [3, [l^ . 

The inclusion of noise of this type may easily destabilize the system, since unrealistically large fluctuations may 
arise even for very small values of a. In the absence of noise, the populations will converge to the stable equilibrium 
given by 



Neq = J (6) 
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(7) 



Figure [2] shows how the populations fluctuate around their average stable values Ngq and Peg, but the fluctuations 
can eventually take the system to the unstable equilibrium N — K and P — 0.0. This negative effect of noise can be 
attenuated by the inclusion of space variables and diffusion or advection, as will be seen in the next section. 

The inclusion of noise does not always have a negative impact on dynamics; in some systems, its effects can 
be particularly important and can lead to unexpected results such as stochastic resonance [l6l . , a phenomenon 
characterized by the enhancement of the response of a system to a periodic driving force in the presence of noise. 
Stochastic resonance has been studied in man y sy stems, including spatially extended systems and pattern-forming 
systems, such as the Swift- Hohenberg equation 1 1811 . Another interesting possibility is the suppression of internal noise 
by the application of an external noise source [l9|. 



IV. THE INCLUSION OF SPACE VARIABLES 



The models considered up to now are in a sense incomplete, since they provide a mean-fleld description of the 
variations of the population in which only global changes are considered. Even though average population densities 
and even some naturally occurring phenomena such as outbreaks can be predicted, many other features of the systems 
are left out, pattern formation being one of these. 



A. Diffusion 

Turing [13] studied how diffusion together with nonlinear local interactions could lead to the emergence of het- 
erogeneities even when starting from homogeneous conditions in a homogeneous environment. Reaction-diffusion 
equations of the form 

^^i{n)+DV^n (8) 

have been used since then to model pattern formation in system on a wide range, from chemical reactions to bacterial 
chemotaxis and animal coat patterns [3, [2l| . Initially, models of plankton spatial heterogeneity and patchiness were 
based on these types of equations, in which the populations are investigated theoretically by some variation on the 
predator-prey models derived from on the Lotka-Volterra model with an added diffusive term. A study of the plankton 



population stability with this model was performed by Steele [22j . Levin and Segel [23| introduced an autocatalytic 
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FIG. 3: Velocity distribution for a 50 x 50 mesh with spacing of 0.25 km, obtained from the seeded eddy model with 12000 
eddies distributed randomly with Rmin ~ 1, Rmax = 16 and a = 0.1 d~^. The arrows indicate the velocity field in each lattice 
point. 

effect in phytoplankton density and differential dispersal rates which favor higher herbivore motility to account for the 
origin of planktonic patchiness. With their hypothesis, they found a transition from a uniform stable state to a new 
steady state in which plant and herbivore are more concentrated in certain regions, depending on the model parameters. 
However, in the case of planktonic populations, diffusion is only important for the movement of plankton on small 
scales (centimeter scales), in which many biological factors are also important. Moreover, the experimental results 
from [3| show that zooplankton is more patchily distributed and indicate that analytical models used for describing 
phytoplankton which are based on a perturbative analysis of the scale-dependent balance between turbulent diffusive 
flux and exponential reproductive growth are inapplicable for explaining the zooplankton heterogeneity over these 
scales. In simulations with models based on reaction-diffusion equations, zoop lankton is found to be less patchily 
distributed than phytoplankton, contrary to experimental observations [23, [2J|- Also, satellite imagery (see Fig. [T]) 
has displayed eddy patterns that represent the phytoplankton biomass and thus demonstrated that plankton patterns 
in the ocean occur on much broader scales and therefore mechanisms other than diffusion should be considered. 



B. Advection 

Experiments performed on large scale in the ocean and satellite imagery have replaced the concept of slow ocean 
currents by the concept of a continuous distribution of more energetic eddies, which may have space scales up to the 
order of 100km. In these larger scales, the dominant form of motions is due to turbulent lateral diffusion generated 
by strong ocean currents and as such is better modelled by an advective term rather than by diffusion, as considered 
in the previous section. However, turbulent motion is computationally expensive and therefore can not be easily 
included in models of population dynamics. This difficulty can be circumvented by the use of a simpler flow model 
such as the two-dimensional non-divergent seeded eddy model [1^ to simulate geostrophic turbulence, which has an 
energy spectrum E{k) cx as proposed by Kraichnan and Charney [l^,i27,]. The stream function ?/j(r) in this model 
is given by 

n 

^(r)=a^(±).i?,^e-(--)V2«^ (9) 

1=1 
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FIG. 4: Energy spectrum of a 250 x 250 with 80000 eddies with 1 < Ri < 64. The fit for larger wavenumbers is in good 
agreement with that expected for geostrophic turbulence {k~^). 



where Ri and are the radius and position of each eddy and a is a calibration constant. The velocity field is then 
given by 



dy ' dx 



(10) 



The probability distribution for a radius size R is taken to be p{R) oc R~^, with lower radius Rmin and higher radius 
Rmax^ E^nd the signs are taken arbitrarily for each eddy. A typical velocity distribution generated by this model is 
seen in Fig.[3l The energy power spectrum, which approximates that expected for geostrophic turbulence is displayed 
in Fig. m (for further discussion, see Sect, ly]). 

Abraham [6] introduced non-diffusive advection to generate patchiness and found that the characteristic spatial 
patterns of phytoplankton and zooplankton are a consequence of the timescales of their response to changes in their 
environment caused by turbulent advection. He extends a model based on logistic growth psj for phytoplankton by in- 
troducing a maturation time for zooplankton and shows that this maturation time is responsible for the determination 
of spatial structure. However, in his model, he does not include grazing saturation and the type of maturation time 
introduced can lead to non-realistic situations, such as growth of zooplankton in the absence of phytoplankton. His 
biological model consists of three coupled differential equations for the carrying capacity (maximum phytoplankton 
concentration attainable), for the phytoplankton and for the zooplankton. The carrying capacity continuously relaxes 
towards a spatially varying background value. He finds that different distributions are due to different response rates 
to changes in the environment caused by turbulent advection. Neufeld et al. [1^ modelled phytoplankton blooms 
which occur both naturally as a result of seasonal changes and as the result of ocean fertilization experiments using 
similar equations, with terms for logistic growth of phytoplankton, grazing by zooplankton and growth and mortality 
for zooplankton. 

In the model we consider [2^, we include both a noise term described in Sect. IIIII and turbulent advection. It 
corresponds to the equations 
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Fn{N,P) ~v - VN, and 
Fp{N, P) - V VP + C(r, t), 



(11) 
(12) 



7 




5^ 

25 
2 

- 1 5 
1 

0^ 





1 

0.9 

□ fi 

□ 7 

□ j6 

□ 5 
0.4 
05 

□ 2 
0.1 
□ 



FIG. 5: Phytoplankton (top) and zooplankton (bottom) distribution. The phytoplankton distribution follows more closely a 
passive tracer distribution (not shown). However, contrary to a passive tracer pattern, the plankton patterns do not remain 
stationary after a transient period. Results from the simulation of Eqs. (|lip and (|12|l on a 250 x 250 lattice corresponding to 
an area of 62.5 x 62.5 km'^. The values of the parameters are a — 1.0 and the others are the same as in Fig. (2] 
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FIG. 6: Transect of the system showing the variations in plankton concentration. Zooplankton (lower curve) has more intense 
variability in small-scales than phytoplankton (upper curve), as observed in experiments. 

where the noise term is characterized by (^(r, t)^(r', i')} = 2[aP{r,t)f5{Y - v')6{t - t'), Fn{N,P) and Fp{N,P) are 
the same as in Eqs. Q and ([5]). We focus on the intrinsic dynamics of the interactions between populations, which can 
be approximated as local in both space and time. Therefore, we assume that the noise is uncorrelated. Colored noise 
in space and time would arise if we considered external perturbations, such as temperature fluctuations. We denote by 
N{r,t) and -P(r,t) the phytoplankton and zooplankton concentrations, respectively. The equations are discretized on 
a two-dimensional mesh with periodic boundary conditions and solved with a fourth order Runge-Kutta method (soj 
for the deterministic phytoplankton growth and with the Milstein scheme for stochastic differential equations for the 
zooplankton growth [31j]. The advective term was solved using the second upwind scheme or donor cell scheme [Hj, a 
method used in problems of computational fluid dynamics. The initial condition is given by homogeneous populations 
of zooplankton and phytoplankton at their stable equilibrium which corresponds to N^g = 1.0 and Peg = 0.225 for 
the parameter values as in Fig. [2l 

In Fig. [HI we display the typical phytoplankton and zooplankton patterns that emerge from the model. The 
patterns presented agree with experimentally observed ones, in which zooplankton display more patchiness than 
phytoplankton [1, H^. It should be pointed out that the patterns generated are not static. In Fig. [51 we present data 
from a transect of the system, which is analogous to what would be observed by a shipborne sampling method. 

The time correlations of the zooplankton is more rapidly decreasing than that of phytoplankton, since the noise 
introduced in the zooplankton growth is delta correlated in time. It is expected that the time correlation of zooplankton 
concentration should decay rapidly to zero. The phytoplankton concentration, on the other hand, is only indirectly 
affected by this uncorrelated noise, through the nonlinear population dynamics and convection. As such, changing 
patterns arise, as can be seen by the fluctuating correlation function. Even though there are fluctuations, the average 
temporal correlation function does not vanish (see Fig. [7]) and oscillates around a flnite value for long times. The 
phytoplankton patterns generally resemble a tracer pattern, with the difference that the latter becomes static in the 
long time limit. Likewise, Fig. [5| displays the spatial correlation functions of the plankton distribution. Again, due 
to the presence of a noise term that is uncorrelated in space, the zooplankton correlations decay faster than the 
phytoplankton correlations. 
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FIG. 7: Time correlation functions for phytoplankton and zooplankton obtained after the system is allowed to relax for a long 
time (r = 2500). After this transient period, the time correlation function of a tracer is equal to 1, which indicates that the 
pattern has become stationary (not shown). C{t) represents the concentration of phytoplankton for the upper curve and of 
zooplankton for the lower curve. 
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FIG. 8: Spatial correlation functions for plankton. C(x) represents the concentration of phytoplankton and of zooplankton. The 
zooplankton correlation function (broken line) decays more rapidly than the phytoplankton correlation function (continuous 
line). Inset: semi- log plot shows that the phytoplankton correlation has an initial exponential decay. 
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FIG. 9: Variance spectra for phytoplankton and zooplankton, obtained from an average over 250 horizontal and 250 vertical 
transects of the system. Phytoplankton displays a steeper slope, as observed in experiments. The exponents obtained are 
—2.7 ± 0.1 for phytoplankton and —1.0 ± 0.1 for zooplankton. 



V. SPECTRAL ANALYSIS 



It is common in studies of plankton distribution to compare the power spectra of their distribution with that of 
physical factors such as temperature. Spectra of sea-surface temperature derived from airborne radiometer measure- 
ments have exponents that vary between —2 and —3. The ones obtained from shipbornc measurements vary widely, 
but are found in a similar range. Gower et al. calculated the power spectrum of phytoplankton concentration 
patterns averaged over all directions and found that the signal variance is related to the inverse wavelength with 
an exponent of —2.92. Initially, it was proposed that this exponent, being close to the —3 exponent expected for 
the energy spectrum of geostrophic turbulence (see Fig. [5]), indicated that phytoplankton was behaving as a passive 
tracer. However, it was pointed out that a conserved passive tracer in geostrophic flow should follow a law [33 |. 
This discrepancy can be understood from the fact that phytoplankton, being microscopic plants, is better described 
as a reactive tracer because they reproduce, grow, die and are eaten. Besides this, the biological time rates are short 
when compared to typical time-scales of quasi-geostrophic eddies. Reactive tracers stirred by geostrophic turbulence 
generally have spectral slopes in the range from —1 to —3. Smith et al. [35|] found phytoplankton spectral slopes of 
about offshore, k^^-^ inshore and k^^ at length scales < 10 km by analyzing satellite images. Abraham @, [sH] 
showed that biological factors such as growth rates have a strong influence on the spectral slope. The ratio of the 
biological and flow timescales is an important factor and fast-growing organisms such as phytoplankton tend to have 
steeper spectral slopes. 

There have been studies that compare the spectra of plankton patchiness to the k^^/^ model of 3D isotropic 



turbulence [37[ . It has been argued in the review by Franks [38| that planktonic data in these studies do not resolve 



the appropriate (small) spatial scales in order to make the comparison. Except in intense turbulence, the small spatial 
scale is characterized by the inertial subrange which lies typically in between scales of about 1 and 100 cm. 

We analyze larger scales at which the water motion is better described by two-dimensional geostrophic turbu- 
lence [1^ [231 • In Fig. [HI we present the plankton spectra obtained from averaging the power spectrum of 500 transects, 
250 in the vertical direction and 250 in the horizontal. We obtain spectral slopes comparable to experimental results, 
with zooplankton having a more white-noise-like spectrum. This is due to the random term in the dynamics that 
leads to a higher variability in zooplankton. The phytoplankton population is also affected by noise, but in an indirect 
way. 
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VI. CONCLUSION 



Even though pattern formation in plankton populations is a long known and well documented phenomena, there 
is still controversy about the mechanisms that give rise to their spatial variability. We believe that simplified models 
capturing the main features of the population dynamics can help to shed light on the dynamics of the actual system. 
The addition of noise allows to have a more realistic description, because zooplankton population is smaller and 
therefore is not as well modelled by a continuous field. Besides this, the interplay between advection and the growth 
dynamics with a random term generate results that agree better with experimental data than models based on 
reaction-diffusion equations. 
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